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ABSTRACT 

We present a detailed investigation of the y-ray emission in the vicinity of the supernova remnant (SNR) W28 
(G6.4— 0.1) observed by the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope. We 
detected significant y-ray emission spatially coincident with TeV sources HESS J1800— 240A, B, and C, located 
outside the radio boundary of the SNR. Their spectra in the 2-100 GeV band are consistent with the extrapolation 
of the power-law spectra of the TeV sources. We also identified a new source of GeV emission, dubbed Source 
W, which lies outside the boundary of TeV sources and coincides with radio emission from the western part of 
W28. All of the GeV y- ray sources overlap with molecular clouds in the velocity range from 0 to 20 km s _1 . 
Under the assumption that the y-ray emission toward HESS J1800— 240A, B, and C comes from jr° decay due 
to the interaction between the molecular clouds and cosmic rays (CRs) escaping from W28, they can be naturally 
explained by a single model in which the CR diffusion coefficient is smaller than the theoretical expectation in the 
interstellar space. The total energy of the CRs escaping from W28 is constrained through the same modeling to be 
larger than ~2 x 10 49 erg. The emission from Source W can also be explained with the same CR escape scenario. 
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1. INTRODUCTION 

Diffusive shock acceleration (DSA) operating at the shocks 
of supernova remnants (SNRs; Reynolds 2008 and references 
therein) is the most likely mechanism to convert the kinetic 
energy released by supernova explosions into high energy parti- 
cles (cosmic rays; CRs) that obey a power-law type distribution. 
Evidence of the CR proton acceleration in SNRs has emerged 
from the detection of GeV y- rays from some SNRs interact- 
ing with molecular clouds such as W51C, W44, and IC 443 
(Abdo et al. 2009b, 2010a, 2010b) by the Large Area Telescope 
(LAT) on board the Fermi Gamma-ray Space Telescope. The 
intense GeV emission from those SNRs is naturally explained 
by 7r° decay produced in inelastic collisions of the acceler- 
ated protons with dense gas. This was recently confirmed by 
the detection of the characteristic spectral feature produced by 
the decay of tt ° s in W44 and IC 443 (Giuliani et al. 2011; 
Ackermann et al. 2013). In DSA theory, CRs accelerated at the 
shock are scattered by self-generated magnetic turbulence. Since 
the highest-energy CRs in the shock precursor at the upstream 
are prone to lack self-generated turbulence, they are expected to 
escape from the shock (Ptuskin & Zirakashvili 2003). However, 
it has been unclear how the CRs escape from SNRs and prop- 
agate into the interstellar medium (ISM) because the interplay 
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among the CRs, the magnetic turbulence, and the surrounding 
environment of SNRs is not well understood. 

If an SNR is in a dense environment, we can expect an 
enhancement of the 7r°-decay y- rays from molecular clouds 
illuminated by the escaping CRs in the vicinity of the SNR 
(Aharonian & Atoyan 1996; Gabici et al. 2007). For example, 
the y-ray emissions near middle-aged SNRs G8.7— 0. 1 and W44 
are naturally explained by the above scenario (Ajello et al. 2012; 
Uchiyama et al. 2012). The energy dependence of the diffusion 
coefficient of the CRs alters their spectrum, which affects the 
spectral shape of the resulting y-ray emission (Aharonian & 
Atoyan 1996; Gabici et al. 2009; Ohira et al. 2011). Thus, we 
can constrain the diffusion coefficient by measuring the wide- 
band y-ray spectrum of the emission around SNRs. 

W28 (also known as G6.4— 0.1) is a mixed-morphology SNR 
whose age is estimated to be (3.5-15) x 10 4 yr (Kaspi et al. 
1993). In this paper, we adopted the same age of 4 x 10 4 yr 
as used in Abdo et al. (2010c). The SNR is located within a 
molecular cloud complex with a mass of 1.4 x 10 6 M e (Reach 
et al. 2005) and interacts with some parts of the cloud, traced 
by the detection of OH (1720 MHz) masers (Frail et al. 1994; 
Claussen et al. 1997, 1999). Observations of molecular lines 
placed W28 at a distance of ~ 1.9 kpc (Velazquez et al. 2002). 
GeV y-ray emission associated with W28 has been detected by 
the LAT and the Gamma-Ray Image Detector (GRID) on board 
AGILE (Tavani et al. 2008). A natural explanation is the decay 
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of jt° s due to the interaction of the cloud and CRs accelerated 
in the SNR (Abdo et al. 2010c; Giuliani et al. 2010). W28 
is considered to have entered the radiative phase (Lozinskaya 
1992) as indicated by optical filaments (Lozinskaya 1974). Thus, 
we can expect that CRs have escaped into the surrounding 
ambient medium. 

H.E.S.S. observations of W28 have revealed four TeV 
y- ray sources (Aharonian et al. 2008): HESS J 1 80 1 — 233, lo- 
cated along the northeastern boundary of W28, and a complex 
of three sources, HESS J1800— 240A, B, and C, located to the 
south, outside the radio boundary. The southern H.E.S.S. sources 
spatially correspond with molecular clouds whose distances are 
consistent with that of W28 (Aharonian et al. 2008), suggesting 
the possibility that their origins are due to runaway CRs from the 
SNR. Thus, this region is one of the best sites to study CR dif- 
fusion. Although Abdo et al. (2010c) detected only one source 
associated with HESS J1800— 240B in the first year of obser- 
vations, there are two LAT sources in the southern region listed 
in the second Fermi-LKT catalog (2FGL; Nolan et al. 2012). If 
the GeV and TeV emissions originate from CRs escaping from 
W28, we can constrain the diffusion coefficient of the particles 
in this region. 

In this paper, we report a detailed analysis of the LAT sources 
surrounding W28, based on four years of data. First, we give 
a brief description of the observations and data selection in 
Section 2. The analysis procedure and results are presented 
in Section 3, along with the spectra of the LAT sources. The 
discussion is given in Section 4 followed by conclusions in 
Section 5. 

2. OBSERVATIONS AND DATA REDUCTION 

The LAT is the main instrument of Fermi, detecting y-rays 
by conversion into electron-positron pairs in the energy range 
from ~20 MeV to >300 GeV (Atwood et al. 2009). It contains a 
high-resolution converter/tracker for direction measurement of 
the incident y-rays, a Csl (TI) crystal calorimeter for energy 
measurement, and an anti-coincidence detector to identify 
the background of charged particles. The LAT has a large 
effective area (~7500 cm 2 on-axis for >2 GeV), a wide field of 
view (~2.4 sr), and a good point-spread function (PSF; the 
68% containment angle at >2 GeV is smaller than ~0?6). 
The on-orbit calibration, event classification, and instrument 
performance are described in detail in Abdo et al. (2009a). 

We have analyzed 4 years of data within the energy range 
2-100 GeV in the vicinity of W28, collected from 2008 August 4 
to 2012 August 18, with a total exposure of ~1.2 x 10 11 cm 2 s 
at 2 GeV. The LAT was operated in sky-survey mode nearly 
continuously. In this observing mode, the LAT scans the whole 
sky every two orbits (~3 hr), obtaining complete sky coverage 
and approximately uniform exposure. 

We used the standard LAT analysis software, ScienceTools 
version v9r30, publicly available from the Fermi Science Sup- 
port Center (FSSC). 14 We used the post-launch instrument re- 
sponse functions (IRFs) P7V6 (Ackermann et al. 2012) and ap- 
plied the following event selection criteria: ( 1 ) events should be 
classified as “Source” class, (2) the reconstructed zenith angles 
of the arrival direction of y- rays should be smaller than 100° to 
minimize the contamination from Earth-limb y- ray emission, 
and (3) only time intervals when the center of the LAT field of 
view is within 52° of the local zenith are accepted, to further 


14 Software and documentation of the Fermi ScienceTools are distributed by 
Fermi Science Support Center at http://fermi.gsfc.nasa.gov/ssc. 
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Figure 1 . Fermi LAT 10-100 GeV count map around W28. The count map is 
smoothed by a Gaussian kernel of a = 0?20. The pixel size is 0?05. The green 
circles and crosses indicate the extended and point sources in the 2FGL catalog 
(Nolan et al. 2012), respectively. The white line from top left to bottom right 
indicates the Galactic plane. 

(A color version of this figure is available in the online journal.) 

reduce the contamination by Earth’s atmospheric emission. The 
y- ray burst GRB 100826A (McEnery & Omodei 2010) occurred 
within the region used for the analysis in this paper. However, 
the event data are not included because they do not satisfy the 
above criteria (3). Thus, we did not need to apply any additional 
time cut. 

Two different tools were used to perform the spatial and 
spectral analysis: gtlike (in binned mode) and pointlike, 
gtlike is a standard maximum likelihood method (Mattox et al. 
1996). pointlike is an alternate binned likelihood technique, 
optimized for characterizing the extension of a source that has 
been extensively tested against gtlike (Lande et al. 2012). 

3. ANALYSIS AND RESULTS 
3.1. Morphological Analysis and Source Position 

Here, we analyzed the morphology of high energy y- ray 
emission in the vicinity of W28. We made a count map in the 
10-100 GeV energy band to take advantage of optimal angular 
resolution and weaker Galactic diffuse emission. Figure 1 shows 
the map in a 10° x 10° region around W28. Emission around 
W28 can be clearly seen against the Galactic diffuse emission. 
There are two LAT sources in the vicinity of W28 in the 2FGL 
catalog (Nolan et al. 2012): 2FGL J1758.8-2402c and 2FGL 
J1800.8— 2400. 

Figure 2 shows a close-up view of the LAT counts map, 
superimposed on the H.E.S.S. significance map (Aharonian et al. 
2008) and the radio image by the Very Large Array (VLA; 
Brogan et al. 2006). Multiple spatial associations are evident, 
allowing GeV and TeV emission to be correlated between the 
northern part of W28 and HESS J 1 801 — 233, between 2FGL 
J1758.8— 2402c and HESS J1800-240C, and between 2FGL 
J 1800.8-2400 and HESS J1800-240B (Aharonian et al. 2008). 
GeV y- rays also overlap with HESS J1800— 240A although 
there is no 2FGL source there. We also found an additional GeV 
source to the west beyond the observed TeV emission. Note 
that this y- ray source cannot be seen clearly below 10 GeV. 
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Figure 2. LAT count map between 10 and 100 GeV around W28 superimposed on (a) H.E.S.S. and (b) VLA contours. The maps have a pixel size of 0?05 and are 
smoothed by a Gaussian kernel of o = 0?20. The inset of each figure shows the simulated LAT PSF with a photon index of 2.5 in the same energy range, adopting 
the same smoothing. Locations of 2FGL catalog sources included in the emission model are indicated with black marks: a circle for W28, a left-hand cross for 2FGL 
J1800.8— 2400, and a right-hand cross for 2FGL J1758.8— 2402c, respectively. A blue x indicates the best-fit position of Source W. Blue diamonds on the left indicate 
Hu regions: G6. 1—0.6 (Kuchar & Clark 1997) and G6. 225— 0.569 (Lockman 1989). The blue diamond on the right is W28A2 (see the text). The white diamond 
indicates the OH maser spot associated with G5. 7 1—0.08 (Hewitt et al. 2009). Green contours in panel (a) show the H.E.S.S. significance map for TeV y-rays at 20, 
40, 60, and 80% levels (Aharonian et al. 2008). Bright TeV spots in the south are HESS J1800— 240A, B, and C as indicated in the figure. Green contours in panel (b) 
indicate the VLA 90 cm image at 5, 10, and 20% of the peak intensity (Brogan et al. 2006). Outer boundaries of SNR W28 and G5.7 1—0.08, as determined from the 
radio images, are drawn as white dashed circles. 

(A color version of this figure is available in the online journal.) 


Here, we will refer to the GeV emissions coincident with the 
two TeV y ray sources as HESS J1800— 240B and 240C, and 
to the western GeV emission as Source W. In the VLA image, 
Source W overlaps with the western shell of W28, whereas 240B 
and 240C overlap with radio emission unrelated to W28: Hii 
region W28A2 (G5. 89-0.39) and SNR G5.71-0.08 (Brogan 
et al. 2006), respectively. 

In order to evaluate the morphology of the GeV emission 
around W28, we fit models to the data using the maximum 
likelihood framework with pointlike. The likelihood is the 
product of the probabilities of observing the y-ray counts 
within each spatial and energy bin for a specified emission 
model. The best parameter values are estimated by maximiz- 
ing the likelihood of the data over the set of models (Mattox 
et al. 1996). The probability density function for the likeli- 
hood analysis includes (1) individual sources detected in the 
2FGL catalog within 15° of W28; (2) the Galactic diffuse emis- 
sion resulting from CR interactions with interstellar medium 
and radiation based on the LAT standard diffuse background 
model, gal_2yearp7v6 .fits available from the FSSC 15 ; and 
(3) an isotropic component to represent extragalactic y- rays 
and residual CR background using a tabulated spectrum written 
in iso_p7v6source.txt, also available from the FSSC. The 
region of interest for the binned maximum likelihood analysis 
based on Poisson statistics 16 was a square region of 16° x 16° in 
Galactic coordinates centered on W28 with a pixel size of 01 1. 

First, we determined the strength of the diffuse y-ray emission 
around W28. To take advantage of the narrower PSF at higher 


15 The model can be downloaded from 

http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html. 

16 As implemented in the publicly available Fermi Science Tools. The 
documentation concerning the analysis tools and the likelihood fitting 
procedure is available from 

http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/. 


energies, we analyzed data from 2 to 100 GeV. To account 
for any effects from nearby sources on the fit, we set free the 
normalization and spectral index of a power-law model applied 
to the Galactic diffuse emission and the locations and spectral 
normalizations of all 2FGL sources within 4° of the direction of 
W28. To represent the emission near Source W more accurately, 
the spectral parameters of W28, HESS J1800— 240B and 240C, 
were also set free. All parameters of the sources beyond 4° were 
fixed to the values in the 2FGL catalog. After the fit, the spectral 
parameters, except for 240B and 240C, were fixed to the values 
obtained in the above analysis. We used a uniform disk model as 
the spatial template of W28, as in the 2FGL catalog. Fitting the 
extension and position of W28 gave consistent results with those 
in the catalog. We substituted the H.E.S.S. significance map of 
HESS J 1 801 — 233 (Aharonian et al. 2008) for the template and 
found that it cannot represent the GeV emission from W28. 

We performed a series of maximum likelihood fits to inves- 
tigate the GeV morphology around W28, adopting a power-law 
spectral form for all sources of interest. First, we added Source 
Was a point source and varied its position. As a result, we ob- 
tained the resulting maximum likelihood value for three point 
sources composed by 240B, 240C, and Source W(L^ p ) with re- 
spect to that for the no-source hypothesis ( Lq ). The likelihood 
ratio, —2 ln(Lo/^ 3 p) (12 degrees of freedom) of ^335, was sub- 
stantially better than that obtained for a model containing only 
the two point sources 240B and C, —2 ln(Lo/L 2 P ) (8 degrees of 
freedom) »293, where Liy, is the likelihood for the two sources 
with the optimized positions (see Table 1). We therefore con- 
cluded that Source W is significantly detected. Note that the 
locations of the other two sources do not significantly differ 
from those in the 2FGL catalog. Then, we added a point source 
at the peak of the TeV emission, 240A. The obtained likeli- 
hood ratio increased by ~29, so we concluded that 240A is 
also significantly detected. In addition, we substituted the four 
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Table 1 

Results of the Morphological Analysis of the y-Ray Emission from the 
Sources in the Vicinity of W28 (2-100 GeV) 


Model 

— 21n (L 0 /Lf 

Additional Degree of Freedom 

W28 

0 

0 

Two point source s b 

293.1 

8 

Three point sources 0 

334.6 

12 

Four point sources d 

363.5 

16 

H.E.S.S. e + Source W 

382.9 

ll f 


Notes. 

a 21n ( Lo/L ), where Lo and L are the maximum likelihood for W28 with the 
2FGL disk template and the additional source component, respectively. 
b 2FGL J1800.8— 2400 and 2FGL J1758.8-2402c in addition to W28 with 
positions free in the optimization. 

c Two point source model plus Source W with optimized position. 
d Three-point-source model plus a point-source model at the peak of FIESS 
J1800-240A. 

e FI.E.S.S. significance map is used as a template for the intensity of the y-ray 
emission. 

f One degree of freedom for the choice of the threshold used for extraction of 
the regions. 

point sources model with a morphological model composed by 
Source W in addition to the H.E.S.S. significance map. The map 
was divided into individual sources 240A, 240B, and 240C with 
separate spectral parameters whose boundaries were determined 
from the apparent TeV morphology. For the latter, we extracted 
the regions above 4 a to avoid the background fluctuations. The 
resulting maximum likelihood ratio of s»383 ( 1 1 degrees of free- 
dom) is larger than that of the four point sources model. Thus, we 
concluded that the best-fit model for the GeV emission is pro- 
vided by including the FI.E.S.S. template and a separate source 
Source W, and used this model for the spectral analysis. Note 
that this result holds even if the extraction threshold for the tem- 
plate extraction from the H.E.S.S. significance map is changed 
by ±lcr. 

Source W is consistent with a newly detected point source. 
The upper limit on its radius, assuming a uniform disk spatial 
model, was 4' at a 68% confidence level. The position of 
Source W in J2000 was obtained as (R.A., deck) = (17 h 58™2, 
— 23°42:3) with an error radius of 05033 at a 68% confidence 
level. We found no other obvious multiwavelength counterparts 
to GeV sources, such as pulsars and blazars, within the positional 
error radius of 0.054 at the 95% confidence region. We tested 
the possibility that Source W is a background active galactic 
nucleus, such as a blazar, which typically has a longer variability 
timescale of the y-ray flux than a few months. Although 
Source W is not bright enough to investigate the variability on 
monthly timescales, we can expect that the test statistics (TS) 
of Source W slowly increase each year if the source is steady. 
Indeed, we found that its TS gradually increases with time. 
However, it remains possible that Source Wmay be a y -ray 
blazar with repetitive yearly activity. We also could not exclude 
the possibility that Source W is a y-ray pulsar due to the lack 
of photon statistics for a pulsation search. 

3.2. Energy Spectrum 

For the spectral analysis of LAT sources in this region, we 
used the maximum likelihood fit tool, gtlike. Each source was 
modeled as a simple power-law function ( dN(E)/dE ) oc E~ a 
for the fit. The obtained spectral index of 2.77 ± 0.06 and 
flux level for W28 were in agreement with the previous result 


Table 2 

Power-law Spectral Indexes and Test Statistics for the LAT Sources Near W28 
in the 2-100 GeV Band 


Name 

Index a 

TS 

HESS J1800-240A 

2.12 ±0.23 ±0.14 

37 

HESS J1800-240B 

2.45 ±0.19 ±0.07 

88 

HESS J1800-240C 

2.38 ±0.23 ±0.17 

51 

Source W 

2.06 ±0.20 ±0.14 

41 


Note. a The first and second uncertainties listed represent the statistical 
and systematic errors, respectively. 


reported by Abdo et al. (2010c) but the spectral energy distri- 
bution is not shown here because we focus on the sources near 
W28. Figure 3 shows the resulting spectral energy distribution 
for each source along with 68% confidence regions. We show 
the best-fit model of HESS J 1 800— 240C obtained by Aharonian 
et al. (2008) as an upper limit for Source W in the TeV range 
(see Figure 3(d)). The obtained spectral indices and TS values 
are shown in Table 2. The spectral index of HESS J1800— 240B 
obtained here is consistent within the uncertainties with Source 
S in Abdo et al. (2010c). 

The LAT spectra of 240B and 240C smoothly connect to 
the H.E.S.S. measurements, while 240A has a slightly harder 
spectral index than the value of 2.55 ±0.18 found by H.E.S.S. 
(Aharonian et al. 2008). Thus, 240A is expected to have a 
spectral break in the GeV to TeV range. However, a simple 
power-law spectrum with a lower flux and softer spectral index 
at the lcr level away from the LAT best-fit values is consistent 
with the H.E.S.S. spectrum. This result is different from that 
presented in Abdo et al. (2010c) where it was shown that 
240 A has a spectral break in the LAT band. This might be 
caused by the difference in source extension or an improved 
understanding of the diffuse background models. Source W also 
can be expected to have a spectral break in the GeV to TeV band 
from Figure 3(d). We fitted the LAT spectrum with a power law 
with an exponential cutoff. However, no clear evidence for a 
break was found from the GeV data alone. 

We considered the systematic errors due to the extraction 
threshold of the H.E.S.S. significance map, the uncertainty 
of the GeV y- ray morphology of W28, the LAT effective 
area, and the modeling of interstellar emission. We evaluated 
the systematic errors associated with the H.E.S.S. map by 
changing the nominal threshold of 4a that we used to extract 
the morphology templates, as explained in Section 3.1, by ±1 <t. 
To account for imperfections in the spatial model describing the 
morphology of W28 (see Figure 2), we divided the uniform disk 
into four quadrants to more accurately represent the morphology 
and estimated systematic errors. The uncertainties in the LAT 
effective area are 5% at 516 MeV and 10% above 10 GeV, 
linearly varying with the logarithm of energy between those 
values (Ackermann et al. 2012). We estimated the systematic 
errors induced in the source flux by repeating the analysis with 
two sets of modified IRFs, where the effective area was scaled 
up and down by its uncertainty. 

In order to evaluate the systematic uncertainties due to the 
interstellar emission model, we compared the results obtained 
using the standard model in Section 3.1 with the results based 
on eight alternative interstellar emission models as performed 
in Ackermann et al. (2013) and de Palma et al. (2013). We 
varied some of the most important parameters of the interstellar 
emission models, namely, the uniform spin temperature used 
to estimate the column densities of interstellar atomic hydrogen 
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Figure 3. GeV-TeV spectral energy distributions of (a) HESS J1800— 240B (b) 240A, (c) 240C, and (d) Source W. The red regions are the 68% confidence range 
of the LAT spectra. The black regions show the combined systematic errors. The black circles show the data points for the H.E.S.S. measurement (Aharonian et al. 
2008). In panel (b), the best-fit model of 240B obtained by Aharonian et al. (2008) is shown as an upper limit in the TeV range. In each panel, the blue solid, dotted, 
and dashed-dotted lines show the model curves with (a) (Z> 27 , <5, r) = (0.5, 0.35, 25 pc), (1, 0.1, 20 pc), and (5, 0.1, 20 pc), respectively; (b) (D 27 , 8, r) = (0.5, 0.35, 
30 pc), (1, 0.1, 21 pc), and (5, 0.1, 21 pc), respectively; (c) ( D 21 , 8 , r) = (0.5, 0.35, 30 pc), (1, 0.1, 24 pc), and (5, 0.1, 24 pc), respectively; and (d) (£> 27 , 8, r ) = (0.5, 
0.35, 25 pc), (1, 0.1, 16 pc), and (5, 0.1, 16 pc), respectively. The gray curve in each panel indicates the upper limit of the y-ray emission produced by the sea of 
Galactic CRs in the same CR-illuminated cloud. 

(A color version of this figure is available in the online journal.) 


( 1 50 K and 1 0 5 K), the vertical height of the CR propagation halo 
(4 kpc and 10 kpc), and the CR source distribution in the Galaxy 
(the pulsar distribution by Lorimer et al. 2006 and the SNR 
distribution by Case & Bhattacharya 1998). We replaced the 
standard isotropic background and Galactic interstellar emission 
models with the alternative ones for the spectral analysis. In this 
procedure, we fixed the isotropic background spectrum while 
the normalization of the interstellar model components were 
left free. The combined systematic errors on the spectral indexes 
and the spectral shapes considering the above uncertainties are 
shown in Table 2 and Figure 3, respectively. 

4. DISCUSSION 
4.1. HESS J1800-240A, B, and C 

Three GeV y-ray sources are found to be spatially coincident 
with the TeV sources HESS J1800-240A, 240B, and 240C. 
They are likely to be steady sources and overlap with molecular 
clouds. This suggests that the GeV-TeV y-ray emission may be 
produced by tt° decay originating from the interaction between 
the clouds and CRs, which were accelerated in and escaped 
from the SNR. Alternatively, the three H.E.S.S. sources may be 
unrelated to W28 but still close to the molecular clouds. Since 
there is no clear TeV counterpart to Source W, we will discuss 
its origin in the following subsection. 


4.1.1. Individual Known Sources 

First, we consider whether the three H.E.S.S. sources may 
be explained by individual sources unrelated to the runaway 
CRs from W28. HESS 1 1800-240A is coincident with two H ii 
regions, G6. 1-0.6 (Kuchar & Clark 1997) and G6.225-0.569 
(Lockman 1989), that contain young massive stars which might 
contribute to the y-ray emission. 240B is associated with 
the clouds containing the ultra-compact Hii region W28A2, 
containing a massive star in a very young phase of evolu- 
tion. W28A2 exhibits energetic bipolar molecular outflows 
(Harvey & Forveille 1988; Acord et al. 1997; Sollins et al. 
2004), which arise from the accretion of matter by the progeni- 
tor star. Abdo et al. (2010c) concluded that GeV y- rays can be 
explained by 7r°-decay and bremsstrahlung. The kinetic energy 
of the outflow is 3.5 x 10 46 erg and the dynamical timescale is a 
few 10 3 yr, with the matter density as high as 10 7 cm~ 3 (Klaassen 
et al. 2006), although there is no model to explain multi-TeV 
particle acceleration in such Hii regions. 240C is spatially co- 
incident with the radio-faint SNR G5. 71— 0.08 (Brogan et al. 
2006), which was suggested to be interacting with a molecular 
cloud due to the detection of an OH maser (Hewitt et al. 2009). 
Its distance is estimated to be either 3.1 or 13.7 kpc based on the 
maser velocity. Therefore, it is not conclusive that the origins of 
the H.E.S.S. sources are stellar objects, or otherwise unrelated 
to SNR W28. 
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4.1.2. Cosmic-Ray Escape Model 

We now explore the possibility that the three H.E.S.S. 
sources are attributed to the 7r°-decay /-rays from molecular 
clouds illuminated by the escaping particles accelerated in W28. 
Attempts to constrain the diffusion coefficient of the runaway 
CRs have been made recently using the GeV-TeV spectrum 
of HESS J1800— 240B and the TeV data with upper limits in 
the GeV band (e.g., Gabici et al. 2010; Li & Chen 2010). We 
can expect that the diffusion coefficient may be more tightly 
constrained by spectral modeling using the GeV-TeV spectra 
of the H.E.S.S. sources presented in this work. 

The CR escape scenario generally assumes that particles 
accelerated in the SNR are gradually released into the ambient 
medium. Here, we assume an energy-dependent release of 
accelerated particles after the time fsT> when the SNR enters 
the Sedov phase (e.g., Gabici et al. 2009; Ohira et al. 2011). 
To estimate the fsT of W28, we adopt the following parameters; 
the typical kinetic energy released by the supernova explosion 
Esn = 10 51 erg and the ejecta mass M e] = 1.4 M 0 . Assuming 
evolution in the uniform intercloud gas with a hydrogen number 
density of 2 cm~ 3 (Reach et al. 2005), the Sedov phase started 
around t = fsx — 310 yr when the radius of W28 was 
r ST - 2.4 pc. 

Let us consider the diffusion process of CRs from SNRs. We 
assume that CRs with a momentum, p, can escape from an SNR 
at a time, t = t esc (p), when the SNR radius reaches R esc (p). 
tesc(p) becomes larger for the CRs with lower momentum and 
is assumed to depend on momentum as a power law, starting 
with Csc(Pmax) = fs'r as follows (Gabici et al. 2009; Ohira et al. 
2011 ): 

tesc(p) = ?st ) • (!) 

\ P max / 

Using the Sedov-Taylor solution and Equation (1), one finds 
R esc (p ) = >'%'[( P / Prmix) 2/5/( . Here, we adopt the maximum 
momentum of the particles p max = 10 15 eV c -1 (reached at r s i ) 
and x — 3, following Gabici et al. (2009) and Ohira et al. (201 1). 
Assuming the age of W28 is 4 x 10 4 yr, the SNR is currently 
releasing CRs with p ~ 0.5 GeV c -1 . 

The momentum spectrum of the runaway CRs integrated over 
the SNR expansion is expected to have the form N esc (p ) oc p~ s 
and s ~ 2, if the maximum momentum of the CRs confined in 
the SNR is a power-law function of time such as Equation (1) 
(e.g., Ptuskin & Zirakashvili 2005). However, the index of 
the spectrum could be different from s = 2, since it depends 
on the time history of acceleration efficiency and maximum 
energy (Ohira et al. 2010; Caprioli et al. 2010). Assuming that 
particles are injected into DSA from the thermal plasma at the 
downstream of the SNR shock (Maikov & Voelk 1995) and 
/ = 3, the value of s is changed to be ~2.2 (Ohira et al. 2010). 
We parameterize the total spectrum of CRs injected into the 
interstellar space as 

N esc (p) = k esc p ~ 22 exp (~p/p max ). (2) 

The spatial distribution of the escaped CRs, n(p, r, t), at a 
time t after the supernova explosion and at a distance r from the 
SNR center is described by the well-known diffusion equation 
for the point source case, and it can be solved using the method 
developed by Atoyan et al. (1995). Considering the escape from 
the surface of the expanding SNR shell, Ohira et al. (2011) 


provided the analytical solution of n(p, r, t) as 
where 

R d (p, t ) = 2^D lsu (p)[t - fesc(p)]- (4) 

T ) ism(t) is the diffusion coefficient of the interstellar medium 
and is often parameterized with a power-law energy dependence 
as below: 

£>,»(,,) =10" D, 7 ( T55 tt_)‘cmV l , (5) 

where Dij is the normalization constant. 

4.1.3. Spectral Modeling to Constrain the Diffusion Coefficient 

To minimize the manual scanning of the parameters for con- 
straining of the diffusion coefficient when modeling the other 
clouds, we first consider the spectrum of the y- ray emission to- 
ward the molecular cloud around HESS J1800— 240B because 
its GeV-TeV spectrum is the best determined among the three 
sources. We adopt 1.9 kpc for the distance from the Earth to 
W28. Only CR protons are taken into account since the leptonic 
emission is unimportant in the case of an electron to proton ratio 
of 0.01, as in the local CR abundance. The y- ray spectrum from 
7T° decays produced by the interaction of protons with ambient 
hydrogen is scaled by a factor of 1 .84 to account for helium and 
heavy nuclei in the target material and CR composition (Mori 
2009). The mass of the cloud Mb responsible for 240B is found 
to be 7.0 x 10 4 M 0 by the NANTEN CO (.1 = 1-0) data for 
the velocity range from 0 to 20 km s 1 (Aharonian et al. 2008). 
Because the distance from the SNR center, r, cannot be well 
determined, we treated it as a free parameter. We adopted the 
center of the radio boundary shown in Figure 2(b) as the SNR 
center. The minimum r is the projected distance 20 pc derived 
from the angular distance between the SNR center and the peak 
of the TeV emission. From the above assumptions, our model 
has four adjustable parameters: r, Dxi, <5, and k esc , which is the 
normalization of the CR spectrum. 

The runaway CR spectrum has two cutoffs whose energies 
are determined as (r — R esc ) 2 /R j and (r + R esc ) 2 /R j from 
Equation (3), and R d depends on D 21 from Equations (4) 
and (5). Then, as the value of Dii becomes small, the cutoff 
energies shift to higher energies. Consequently, for a small value 
of D27, the model spectrum of the y- ray emission conflicts with 
the observed GeV spectrum, and the lower limit obtained is 
~0.5. In this case, r must be almost the same value as the 
projected distance, the minimum r, because the cutoff energies 
also move to higher energies with increases in r. On the other 
hand, 8 is constrained to be 0.20-0.35 to fit the spectrum above 
the higher cutoff energy as shown in Table 3. Figure 3(a) shows 
the y- ray model curve with D27 = 0.5, 8 = 0.35, and r — 25 pc. 
The y- ray emission from the Galactic CRs, calculated based 
on the proton spectrum in Dermer (1986), is also plotted to 
show the expected background emission. Even if the interstellar 
emission model does not reproduce the background emission 
from local clouds completely, it is not a large fraction of 
the observed residual emission. The upper limit on D 21 is 
constrained by the amount of escaped CRs, W v , because the 
normalization of the spectrum of the particles depends on 
kesc/Rd oc k esc /Dl^ from Equations (2), (3), (4), and (5). Here, 
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Table 3 

The Values of the Diffusion Coefficient of the Escaping CRs and the Distance from the SNR Center to Each Molecular Cloud Obtained 
by the Spectral Modeling for HESS J1800— 240A, B, C, and Source W 


Z>27 

S 

240A 

240B 

240C 

Source W 

r (pc) 

W„ (10 49 erg) 

r (P c ) 

W p (10 49 erg) 

r (pc) 

W p (10 49 erg) 

r (pc) 

W p (10 49 erg) 

0.5 

0.35 

30^40 

3.3-10 

25 

1.0 

30 

4.3 

25 

10 


0.3 

25-35 

2. 3-4.9 

20-25 

0.5-1. 0 

24-30 

2. 1-4.8 

20-25 

5.4-10 


0.25 

21-35 

1.3-4. 2 

20-25 

0.5-0.8 

24 

2.1 

16-25 

3 .2-7.5 


0.2 

21-30 

1. 0-2.3 

20 

0.4 

24 

1.9 

16-25 

3.2-54 


0.15 

21-30 

0.8-2.0 

N/A 

N/A 

N/A 

N/A 

16-20 

2.6-3.2 


0.1 

21-25 

0 

1 

b 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

1 

0.35 

40-45 

8.2-10 

30-35 

2. 1-2.7 

35 

8.1 

N/A 

N/A 


0.3 

35-45 

5.9-10 

20-35 

1. 1-2.5 

24-40 

4.3-10 

N/A 

N/A 


0.25 

30^45 

3.9-8.2 

20-30 

0.9-1. 6 

24-35 

3.2-64 

25 

10 


0.2 

25^45 

2. 3-7.5 

20-30 

0.8-1. 5 

24-35 

2.1-6.1 

20-25 

7.5-10 


0.15 

21-40 

1 .6-4.7 

20-25 

0.6-0.9 

24-30 

2. 1-3.8 

16-25 

54-7.5 


0.1 

21-40 

1. 1-4.9 

20 

0.5 

24 

2.0 

16 

3.9 

5 

0.35 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 


0.3 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 


0.25 

N/A 

N/A 

30-60 

7.5-10 

N/A 

N/A 

N/A 

N/A 


0.2 

N/A 

N/A 

20-60 

5.4-10 

24 

10 

N/A 

N/A 


0.15 

30-45 

6.4 

20-60 

4.0-10 

24-35 

10 

N/A 

N/A 


0.1 

21-50 

4.3-64 

20-45 

3.2-54 

24-35 

8.1 

N/A 

N/A 


Note. The spectra of all sources can be reproduced with a single diffusion coefficient with D 21 — 0.5 and S — 0.35. 


we assume the upper limit of W p is 10 50 erg, which is 10% of 
£sn- We try to fit the y-ray spectrum with ££7 = 1, 5, and 
10, respectively. ££7 is found to be below ~5 (see Figure 3(a)) 
since the required W p exceeds the upper limit using D 27 = 10. 
The value of S is constrained to be slightly smaller (0.1-0.25) 
than with ££7 = 0.5. A similar result was obtained for IC 443 by 
Torres etal. (2010). The variation of <5 for0.5 < ££7 < 5 is given 
in Table 3. The lower limit of W p is obtained with ££7 = 0.5 and 
the minimum r, as the spectral normalization is proportional to 
k esc /D 2 y r from Equations (2)-(5). Using S — 0.20, the resulting 
value of W p is 0.4 x 10 49 (M B /7.0 x 10 4 M 0 ) _i erg. 

To investigate whether the emission from 240A and 240C 
can be interpreted within the same CR escape scenario, we 
modeled the y -ray spectra considering the range in the diffusion 
coefficient obtained for 240B. The minimum r for 240 A and 
240C are derived from the projected distances of 2 1 pc and 24 pc, 
respectively. Their spectra can be fitted with models using the 
same diffusion coefficient, except for the case of £>27 = 5 (see 
Table 3). The model curves with ££7 = 0.5 and 8 — 0.35 for 240A 
and 240C are shown in Figures 3(b) and (c), respectively, with 
the allowed ranges of r for each combination of the diffusion 
coefficient and are shown in Table 3. The masses of the clouds 
Ma and M c toward 240A and 240C are estimated to be about 
2.3 x 10 4 M o and 1.4 x 10 4 M 0 using NANTEN CO ( J = 
1-0) data by Aharonian et al. (2008) and Nicholas et al. (2012), 
respectively. Using ££7 = 0.5 and 8 = 0.20, we obtained the 
lower limits of W p for 240 A and 240C as 1.0 x 10 49 (Ma/ 
2.3x 10 4 Mo)” 1 erg and 1.9 x 10 49 (M c / 1.4 x 10 4 M o )- 1 erg, 
respectively. The different W p resulting from assuming the same 
diffusion coefficient among the three H.E.S.S. sources may be 
explained by different fractions of the cloud mass responsible 
for the y-ray emissions. With ££ 7 = 5, 8 is constrained to 
be 0.1-0.15. The model curve with <5 = 0.1 is also shown in 
Figures 3(b) and (c). 

By combining the results for the three H.E.S.S. sources, 
the diffusion coefficient for the runaway CRs from W28 is 


constrained to be 0.5 < £>27 < 5 and 0.1 < 8 < 0.35, with 
a negative correlation between them as shown in Table 3. The 
lower limit of W p is obtained to be 1 .9 x 10 49 erg. The value of r 
becomes larger as either £>27 or 8 increase. However, all clouds 
cannot be located beyond three times their projected distances 
(see Table 3). 

4.1.4. Uncertainty of the Diffusion Coefficient 

We also evaluate the uncertainty of the diffusion coefficient 
by changing the values of the following parameters; the age of 
W28, the explosion energy, and the time when the Sedov phase 
starts. Doubling the SNR age changes the range of D 21 to be 
0.5-10 while halving the SNR age changes it to be 0.1-1. To 
take into account the uncertainty in £sn> we adopt 4 x 10 5 ° erg 
as the lower limit because this value was obtained by Rho & 
Borkowski (2002) with the caution of possible underestimation. 
Assuming that the upper limit of W p is 10% of the £sn, £>27 is 
constrained to be less than 1. In this case, 8 slightly changes 
to the harder value of 0. 1-0.2 at D 21 = 1 than those obtained 
above (see Table 3). The ambient matter density may also be 
different from 2 cm~ 3 . If W28 originates from a core-collapse 
supernova, the progenitor makes a cavity into which the SNR 
initially expands, and the time f ST when CRs start to escape will 
change. Adopting a matter density of 0.1 cm~ 3 changes tsr to 
~840 yr. Although the spectral break shifts to higher energy, the 
model curve with each value of ££7 still does not conflict with 
the observed spectrum. In summary, the diffusion coefficient 
does not depend strongly on the choice of the parameters. 

Moreover, there is uncertainty of the ejecta mass. If we 
adopt M e j = 7 Mq, which is likely in core-collapse supernovae, 
tsT becomes ~ 1 100 yr which is longer than that with M e j = 
1.4 M q . Consequently, the y- rays around 2 GeV could be no 
longer explained by the CR escape scenario because, based on 
Equation (1), particles with the momentum lower than about 
10 GeV cannot escape. Thus, it is possible to reproduce all 
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components of the y- ray emission if the ejecta mass is small but 
it would be difficult with masses larger than 7 M G . 

4.2. W28 

The y-ray emission from the northeastern part of W28 may 
also originate from the CR escape scenario. We tried to model 
the y-ray spectrum using the same scenario that we applied to 
theH.E.S.S. sources. Here, r was fixed to the SNR radius of 13 pc 
from the radio boundary. The mass of the cloud responsible for 
the emission Mne is estimated to be 5 x 10 4 M G by Aharonian 
et al. (2008). The model curves with £>27 = 0.5 and 1 reasonably 
reproduce the observed spectrum when 8 = 0.35. Thus, the 
emission is consistent with the CR escape scenario. However, 
there are also two other possible scenarios to explain the y- rays: 
(1) CRs remaining at the SNR shell contribute to the emission 
because this place is the interaction site of W28 with molecular 
clouds, and (2) the emission comes from the reaccelerated CRs 
in the crushed cloud (Uchiyama et al. 2010). In the first scenario, 
y- rays above 2 GeV are emitted by particles above 10 GeV c -1 , 
which should have escaped from the SNR shell because R csc for 
these particles is ~1 1 pc, which is smaller than the SNR radius. 
However, R esc is proportional to M e • and n {) , and given 

the uncertainties in these quantities, may be larger than 13 pc. 
Therefore, we cannot rule out that these CRs remain within the 
SNR and do not use the model results for the northeastern part 
of W28 to further constrain the diffusion coefficient. 

4.3. Source W 

Source W is separated from the TeV emission, having no 
obvious counterpart except for the radio emission of W28. One 
possibility is that the y-rays come from the particles accelerated 
and confined in the SNR shell. The different y- ray spectral 
shape compared with the northeastern boundary of W28, shown 
in Section 3.2, might be due to the difference of the particle 
spectrum if the y-rays from the northeastern region also come 
from the CRs confined in the SNR shell. Since the environment 
around Source W is expected to be more tenuous than that of the 
northeastern boundary (see Figure 2 in Aharonian et al. 2008), 
the effects of damping of magnetohydrodynamic turbulence 
due to ion-neutral collisions (Drury et al. 1996) could operate 
differently, causing a break in the spectrum of accelerated 
particles at higher energy (>100 GeV). A harder radio spectral 
index is found around Source W (Dubner et al. 2000), which 
may support this scenario. 

Another possibility is that the emission originates from the 
interaction between molecular clouds and the runaway CRs from 
W28. Source W overlaps clouds in the velocity range from 0 
to 20 km s' 1 (see Figure 2 in Aharonian et al. 2008). The total 
cloud mass, Mw, is estimated to be about 3500 M e from the 
NANTEN CO (7= 1-0) data (Mizuno & Fukui 2004; Takeuchi 
et al. 20 1 0), assuming the same velocity range and a cloud radius 
of ~3 pc derived from the upper limit of the extension of the 
y- ray emission. We perform modeling of the y- ray spectrum 
in the same manner as in Section 4.1.3, considering the range 
of the diffusion coefficient obtained there. The lower limit of r is 
the projected distance of 16 pc. Using D^i = 5, the model cannot 
reproduce the spectrum, as shown in Figure 3(d) and Table 3, 
because the required W p exceeds 1 0 50 erg. Therefore, D27 is 
constrained to be the smaller range 0.5-1 than that obtained for 
the three H.E.S.S. sources. On the other hand, <5 is also tightly 
constrained to be 0.1-0.25 with /A? = 1. The lower limit of W p 
is 3.2 x 10 49 (M w /3.5 x 10 3 M G ) _1 erg, with £>27 = 0.5 and 


S — 0.2. From these results, the four y- ray sources around W28 
can all be explained by CR escape from the SNR with a single 
diffusion coefficient with £>27 — 0.5 and S ~ 0.35. 

4.4. In terpretation of the Obtained Diffusion Coefficient 

The obtained values for £>27 and 8 are smaller than those 
based on the Galactic CR propagation model: I ) 21 ~ 10 and 
<5 ~ 0.6 (Ptuskin et al. 2006; Delahaye et al. 2008). £>27 is also 
smaller than that of W44 derived by Uchiyama et al. (2012), 
where 1 < £>27 < 30. For an SNR in a dense environment, 
£>27 is expected to be as low as ~0.1 (Ormes et al. 1988) and 
to depend on the magnetic field strength in the molecular cloud 
into which the runaway CRs propagate (Gabici et al. 2007). £>27 
and 8 would also be affected by the amplification of Alfven 
waves generated by the escaping CRs (Fujita et al. 2011). The 
value of 8 depends on the assumption of the spectral index of 
the accelerated particles in the SNR, which is related to the 
time history of acceleration efficiency and maximum energy, 
as mentioned above. The hard spectrum of the runaway CRs 
might indicate a harder spectral index than 2.0 for the CRs 
accelerated in W28, and it is suggested to be 1.7 based on the 
radio synchrotron spectrum (Dubner et al. 2000). Nonlinear 
shock modification caused by the efficient CR acceleration can 
produce such a hard spectrum (e.g.. Baring et al. 1999). 

The value of 8 is also strongly dependent on assumptions 
about the evolution of the accelerator. If an SNR continuously 
accelerates particles even after the free expansion phase, the 
spectral index of the escaped particles above the cutoff energy 
s' is .v+<5. This is smaller by a factor of 8/2 than in the case 
of an impulsive source (Aharonian et al. 2008), which is very 
similar to our assumption for W28, i.e., the value of 8 becomes 
larger by a factor of 1.5 than our result. However, it is difficult 
to adapt this scenario to W28. The SNR is considered to 
be in the radiative phase with a shock velocity of less than 
100 km s -1 (Rho & Borkowski 2002) and cannot produce multi- 
TeV particles. Therefore, the value of 8 around W28 must be 
smaller than that of the Galactic CR propagation model. 

5. CONCFUSION 

We analyzed the GeV y- ray emission in the vicinity of SNR 
W28 using fours years of FAT data. We detected GeV y-rays 
spatially coincident with the TeV sources HESS J1800— 240A, 
B, and C, located south of the radio boundary of W28. Their 
spectra in the 2-100 GeV band are consistent with the extrapo- 
lation of power-law emission from the TeV y- ray sources. We 
also detected GeV emission from Source W, located outside 
the boundary of the TeV emission and coinciding with radio 
emission from the western shell of W28. All of the GeV y- ray 
sources overlap with the molecular clouds in the velocity range 
from 0 to 20 km s -1 . 

Assuming that the y- ray emissions from the three H.E.S.S. 
sources are due to the decay of 7r°s produced by the interac- 
tion of the molecular clouds with CRs escaping from W28, 
the GeV-TeV spectra can naturally be explained by a single 
model. We constrain the diffusion constant at 10 GeV c -1 and 
the power-law index of the energy dependence to be 0.5-5 x 
10 27 cm 2 s' 1 and 0.1-0.35, respectively, with a negative corre- 
lation between them. These values are smaller and harder than 
those of the Galactic CR propagation model. Considering the 
masses of the molecular clouds responsible for the emission, 
the lower limit on the total energy of the escaped CRs is con- 
strained to be ~2 x 10 49 erg, in agreement with the conjecture 
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that SNRs are the main sources of the Galactic CRs. The y-rays 
from Source Wean be also interpreted to be the emission orig- 
inating from the interaction of the runaway CRs and molecular 
clouds with the same diffusion coefficient as obtained for the 
H.E.S.S. sources. 

The research of D.F.T. and G.R has been done in the 
framework of the grant AYA20 1 2-39303 and iLINK20 1 1-0303 . 
D.F.T. was additionally supported by a Friedrich Wilhelm Bessel 
Award of the Alexander von Humboldt Foundation. 

The Fermi LAT Collaboration acknowledges generous ongo- 
ing support from a number of agencies and institutes that have 
supported both the development and the operation of the LAT as 
well as scientific data analysis. These include the National Aero- 
nautics and Space Administration and the Department of Energy 
in the United States; the Commissariat a l’Energie Atomique and 
the Centre National de la Recherche Scientifique/Institut Na- 
tional de Physique Nucleaire et de Physique des Particules in 
France; the Agenzia Spaziale Italiana and the Istituto Nazionale 
di Fisica Nucleare in Italy; the Ministry of Education, Culture, 
Sports, Science and Technology (MEXT), High Energy Accel- 
erator Research Organization (KEK), and Japan Aerospace Ex- 
ploration Agency (JAXA) in Japan; and the K. A. Wallenberg 
Foundation, the Swedish Research Council, and the Swedish 
National Space Board in Sweden. 

Additional support for science analysis during the operations 
phase is gratefully acknowledged from the Istituto Nazionale di 
Astrofisica in Italy and the Centre National d’Etudes Spatiales 
in France. 

REFERENCES 

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009a, APh, 32, 193 
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009b, ApJL, 706, LI 
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, Sci, 327, 1 103 
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010b, ApJ, 712, 459 
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010c, ApJ, 718, 348 
Ackermann, M„ Ajello, M„ Allafort, A., et al. 2012, APh, 35, 346 
Ackermann, M„ Ajello, M„ Allafort, A., et al. 2013, Sci, 339, 807 
Acord, J. M„ Walmsley, C. M„ & Churchwell, E. 1997, ApJ, 475, 693 
Aharonian, F. A., Akhperjanian, A. G., Bazer-Bachi, A. R.. et al. 2008, A&A, 
481,401 

Aharonian, F. A., & Atoyan, A. M. 1996, A&A, 309, 917 
Ajello, M„ Allafort, A., Baldini, L„ et al. 2012, ApJ, 744, 80 
Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071 
Atoyan, A. M„ Aharonian, F. A., & Volk, H. J. 1995, PhRvD, 52, 3265 
Baring, M. G., Ellison, D. C., Reynolds, S. R, Grenier, I. A., & Goret, P. 
1999, ApJ, 513, 311 

Brogan, C. L., Gelfand, J. D., Gaensler, B. M., Kassim, N. E., & Lazio, T. J. W. 
2006, ApJL, 639, L25 

Caprioli, D., Amato, E., & Blasi, P. 2010, APh, 33, 160 
Case, G. L., & Bhattacharya, D. 1998, ApJ, 504, 761 

Claussen, M. J., Frail, D. A., Goss, W. M., & Gaume, R. A. 1997, ApJ, 
489, 143 

Claussen, M. J., Goss, W. M., Frail, D. A., & Desai, K. 1999, ApJ, 
522, 349 


de Palma, F.. Brandt, T. J., Johannesson, G., Tibaldo, L., & Fermi LAT 
Collaboration, 2013, in Proc. of 4th International Fermi Symposium Proc., 
ed. N. Omodei, G. Senatore, T. Brandt, & C. Wilson-Hodge (Stanford, CA: 
Stanford Univ.), 1 72 

Delahaye, T., Lineros, R., Donato, F., Fornengo, N., & Salati, P. 2008, PhRvD, 
77, 063527 

Dermer, C. D. 1986, A&A, 157, 223 

Drury, L. O’C., Duffy, R, & Kirk, J. G. 1996, A&A, 309, 1002 
Dubner, G. M., Velazquez, P. F., Goss, W. M., & Holdaway, M. A. 2000, AJ, 
120, 1933 

Frail, D. A., Goss, W. M„ & Slysh, V. I. 1994, ApJL, 424, LI 1 1 
Fujita, Y., Takahara, F„ Ohira, Y., & Iwasaki, K. 2011, MNRAS, 
415, 3434 

Gabici, S., Aharonian, F. A., & Blasi, P. 2007, Ap&SS, 309, 365 
Gabici, S., Aharonian, F. A., & Casanova, S. 2009, MNRAS, 396, 1629 
Gabici, S., Casanova, S., Aharonian, F. A., & Rowell, G. 2010, in SF2A-2010: 
Proceedings of the Annual Meeting of the French Society of Astronomy 
and Astrophysics, ed. S. Boissier, M. Heydari-Malayeri, R. Samadi, & D. 
Valls-Gabaud (Paris: French Society of Astronomy and Astrophysics), 313 
Giuliani, A., Cardillo, M., Tavani, M., et al. 2011, ApJL, 742, L30 
Giuliani, A., Tavani, M., Bulgarelli. A., et al. 2010, A&A, 516, LI 1 
Harvey, P. M„ & Forveille, T. 1988, A&A, 197, L19 
Hewitt, J. W., Yusef-Zadeh, F„ & Wardle, M. 2009, ApJL, 706, L270 
Kaspi, V. M., Lyne, A. G., Manchester, R. N., et al. 1993, ApJL, 409, L57 
Klaassen, P. D., Plume, R., Ouyed, R., von Benda-Beckmann, A. M., & Di 
Francesco, J. 2006, ApJ, 648, 1079 
Kuchar, T. A., & Clark, F. O. 1997, ApJ, 488, 224 
Lande, J., Ackermann, M., Allafort, A., et al. 2012, ApJ, 756, 5 
Li, H.. & Chen, Y. 2010, MNRAS, 409, L35 
Lockman, F. J. 1989, ApJS, 71, 469 

Lorimer, D. R.. Faulkner, A. J., Lyne, A. G., et al. 2006, MNRAS, 372, 777 
Lozinskaya, T. A. 1974, SvA, 17, 603 

Lozinskaya, T. A. 1992, Supernovae and Stellar Wind in the Interstellar Medium 
(Melville, NY: AIP) 

Maikov, M. A., & Voelk, H. J. 1995, A&A, 300, 605 

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396 

McEnery, J., & Omodei, N. 2010, GCN, 11155 

Mizuno, A., & Fukui, Y. 2004, in ASP Conf. Ser. 317, Milky Way Surveys: 
The Structure and Evolution of our Galaxy, ed. D. Clemens, R. Shah, & T. 
Brainerd (San Francisco, CA: ASP), 59 
Mori, M. 2009, APh, 31, 341 

Nicholas, B. P„ Rowell, G., Burton, M. G„ et al. 2012, MNRAS, 419, 251 
Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012, ApJS, 199, 31 
Ohira, Y„ Murase, K„ & Yamazaki, R. 2010, A&A, 513, A17 
Ohira, Y„ Murase, K„ & Yamazaki, R. 201 1 . MNRAS, 410, 1577 
Ormes, J. F„ Ozel, M. E„ & Morris, D. J. 1988, ApJ, 334, 722 
Ptuskin, V. S., & Zirakashvili, V. N. 2003, A&A, 403, 1 
Ptuskin, V. S., & Zirakashvili, V. N. 2005, A&A, 429, 755 
Ptuskin, V. S.. Moskalenko, I. V., Jones. F. C., Strong, A. W., & Zirakashvili, 
V. N. 2006, ApJ, 642, 902 

Reach, W. T„ Rho, J., & Jarrett, T. H. 2005, ApJ, 618, 297 

Reynolds, S. P. 2008, ARA&A, 46, 89 

Rho, J., & Borkowski, K. J. 2002, ApJ, 575, 201 

Sollins, P. K„ Hunter, T. R„ Battat, J., et al. 2004, ApJL, 616, L35 

Takeuchi, T„ et al. 2010, PASJ, 62, 557 

Tavani, M., Barbiellini, G., Argan, A., et al. 2008, NIMPA, 588, 52 
Torres, D. F., Marrero, A. Y. R., & de Cea Del Pozo, E. 2010, MNRAS, 
408, 1257 

Uchiyama, Y., Blandford, R. D., Funk, S., Tajima, H., & Tanaka, T. 2010, ApJL, 
723, L122 

Uchiyama, Y., Funk, S., Katagiri, H., et al. 2012, ApJL, 749, L35 
Velazquez, P. F., Dubner, G. M., Goss, W. M., & Green, A. J. 2002, AJ, 
124,2145 


9 


